**---------------- estimate production function for each 2-digit industry

******* create variables
** calculate market share, output prices. 
gen output_nominal=V207
bysort cic_adj year: egen ind_output_nominal=sum(output_nominal)
gen ms=output_nominal/ind_output_nominal

gen price=output_nominal/amount
gen p=ln(price)

** log quantity output
gen y=ln(amount) 


** adjusted material input
merge m:1 cic_adj using "${deflators_data}\benchmark_input_deflator_rv.dta"
drop _merge 
gen material_nominal=F356-F340   //intermediate inputs - financial cost
gen material_real=.
forvalues i = 1998(1)2007 {
replace material_real=(material_nominal/InputDefl`i'_updated)*100 if year==`i'
}
gen m=ln(material_real)



**input share
gen alpha_m= material_nominal/output_nominal


** adjusted capital stock: zhangyifan's method
sort id year
merge 1:1 id year using "${processed_data}\real_capital.dta"
drop if _merge==2
drop _merge
gen k=ln(capital_real)



**employment
gen l=log(V210)

**trade
gen e=trade_custom


***generate interation terms
gen p2=p^2 
gen ms2=ms^2
gen pms=p*ms 
gen l2=l^2 
gen m2=m^2 
gen k2=k^2
gen lm=l*m
gen lk=l*k 
gen km=k*m 
gen lmk=l*m*k 
gen l_lag=L.l
gen k_lag=L.k 
gen m_lag=L.m 
gen p_lag=L.p 
gen e_lag=L.e 
gen ms_lag=L.ms  
gen p2_lag=L.p2  
gen ms2_lag=L.ms2 
gen p_lagms_lag=p_lag*ms_lag


drop if y==. | l==. | m==. | e==. | p==. | ms==. 
gen region=substr(B056,1,4)
xi i.prod_id, pre(dma)
xi i.region, pre(dmb)
xi i.year, pre(dmc)

* OLS REGRESSION FOR STARTING VALUES
reg y  l m k  p ms e p2 ms2 pms  dmaprod* dmbregion* dmcyear*
gen bconsols=_b[_cons]
gen blols=_b[l]
gen bmols=_b[m]
gen bkols=_b[k]

gen bpols=_b[p]
gen bmsols=_b[ms]
gen beols=_b[e]
gen bp2ols=_b[p2]
gen bms2ols=_b[ms2]
gen bpmsols=_b[pms]

*-------------------------------------------------------------------
*step 1
foreach var of varlist l m k l2 m2 k2 lm lk km lmk p ms p2 ms2  {
gen rege_`var'=e*`var'
}
xi: reg y l m k l2 m2 k2 lm lk km lmk  rege_* p ms e p2  ms2  pms  dmaprod* dmbregion* dmcyear*  
predict phi
gen phi_lag= L.phi
drop dm*
gen const=1
sort fid year
drop if l_lag==.
drop if m_lag==.
drop if phi==.
drop if phi_lag==.



*----------------MATA PROGRAM------------------*
clear mata
mata:
void GMM_linear(todo,betas,crit,g,H)
{
	PHI=st_data(.,("phi"))
	PHI_LAG=st_data(.,("phi_lag"))
	Z=st_data(.,("const","l_lag","m_lag","k","p_lag","ms_lag", "e_lag","p2_lag","ms2_lag","p_lagms_lag"))
	X=st_data(.,("const","l","m","k", "p", "ms", "e","p2", "ms2","pms"))
	X_lag=st_data(.,("const","l_lag","m_lag","k_lag", "p_lag", "ms_lag", "e_lag","p2_lag","ms2_lag","p_lagms_lag"))
	Y=st_data(.,("y"))
	E_lag=st_data(.,("e_lag"))
	RC_lag=st_data(.,("RC_L"))
	C=st_data(.,("const"))

       OMEGA=PHI-X*betas'
       OMEGA_lag=PHI_LAG-X_lag*betas'        // first-order terms
       OMEGA_lag2=OMEGA_lag:*OMEGA_lag       // second-order terms
       OMEGA_lag3=OMEGA_lag2:*OMEGA_lag      // third-order terms
       OMEGA_lag_pol=(C,E_lag,RC_lag,OMEGA_lag) 
  	g_b = invsym(OMEGA_lag_pol'OMEGA_lag_pol)*OMEGA_lag_pol'OMEGA
	XI=OMEGA-OMEGA_lag_pol*g_b
	crit=(Z'XI)'(Z'XI)
}

void MODEL_linear()
	{
	initialvalue=st_data(1,("bconsols","blols","bmols","bkols","bpols","bmsols","beols","bp2ols","bms2ols","bpmsols"))
S=optimize_init()
optimize_init_evaluator(S, &GMM_linear())
optimize_init_evaluatortype(S,"d0")
optimize_init_technique(S, "nm")
optimize_init_nmsimplexdeltas(S, 0.1)
optimize_init_which(S,"min")
optimize_init_params(S,initialvalue)
p=optimize(S)
p
st_matrix("beta",p)
}
end



mata MODEL_linear()
gen betaconst=beta[1,1]
gen betal=beta[1,2]
gen betam=beta[1,3]
gen betak=beta[1,4]
gen betap  =beta[1,5]
gen betams =beta[1,6]
gen betae  =beta[1,7]
gen betap2  =beta[1,8]
gen betams2  =beta[1,9]
gen betapms  =beta[1,10]
contract cic2 betaconst- betapms 
drop _freq

